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Abstract. We construct, on continuous Qji finite elements over Cartesian 
meshes, an interpolation operator that does not increase the total variation. 
The operator is stable in L 1 and exhibits second order approximation prop- 
erties. With the help of it we provide improved error estimates for discrete 
minimizers of the total variation denoising problem and for total variation 
flows. 



1. Introduction 

The approximation of weakly differentiable functions by polynomials is a very 
useful tool when trying to understand the behavior of such functions. In fact, 
S.L. Sobolev himself [25] used a kind of averaged Taylor polynomial to discuss 
equivalent norms in the celebrated spaces that now bear his name. Polynomial 
approximation also plays a crucial role in numerics, as it is the basis of the anal- 
ysis of finite element methods and discretization techniques of partial differential 
equations. It is of no surprise then, that there exists several different constructions 
and that the approximation properties of such polynomials are well studied. We 
refer the reader to El [101 EEU H2 EB El] for a by no means exhaustive list that 
refers, mostly, to such results in the finite element method context. 

The quasitotality of the aforementioned works construct an operator that is well 
defined in L 1 (J7) (see Section [5] for notation) and is stable in L p (il) for p > 1. In 
addition, their approximation properties are studied when the functions belong to 
certain Sobolev spaces. However, in problems dealing with free discontinuities or 
connected with geometric features such as minimization of coarea; mean curvature 
flows and others, it is most natural to work in the space BV(Q) and so the need 
for an approximation theory by piecewise polynomials in this space arises. 

In prior work [4] we dealt with such issue by first regularizing the function and 
then taking the Lagrange interpolant of the regularization. By properly choosing 
the regularization parameter as a function of the mesh size, we were able to ob- 
tain the properties that served our purposes. However, much simpler arguments 
and better results would have been possible if a TV diminishing interpolant was 
available. In other words, we wonder about an operator 

(TVD) U h :L 1 (Q)^X h : [ \VU h w\ < [ \Vw\, V» e By(Jl), 

Jn Jn 
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where we denote by Xt the finite clement space; note that the stability constant in 
(|TVDj) is 1. The purpose of this work is to construct such an operator under the, 
rather stringent, geometric assumption that the underlying mesh is Cartesian and 
the finite element space is made of continuous piecewise multilinear elements. 

We organize our presentation as follows. Notation, functional and discrete spaces 
are introduced in Section [2j The core of our discussion is Section |3j where we 
construct our TV diminishing interpolant and prove its main properties for periodic 
functions. The extension to non-periodic functions is in Section 2J Finally, some 
applications are discussed in Section [5] namely we provide error estimates for total 
variation minimization and total variation flows. 

2. Notation and preliminaries 

In what follows we consider either = S 1 x S 1 (the two dimensional torus) or 
f2 = (0, l) 2 (the unit square). Our arguments can be extended without difficulty to 
more space dimensions. If z £ R 2 we denote its coordinates by z = (z , z 2 ). 

We denote by L p (tt) with p £ [1, oo] the space of Lebesgue integrable functions 
with exponent p. Sobolev spaces will be denoted by Wp(ft), where s € R is the 
differentiability order. Whenever X is a normed space, we denote by || • \\x its 
norm. 

We say that a function w £ BV(il) (is of bounded variation), if w £ i 1 (51) and 
its derivative in the sense of distributions is a Radon measure. To simplify notation, 
we will use Ww to denote such measure. We define the total variation (TV) by 

(2.1) TVH = f |Vtu|, 

Jn 

and we endow the space BV(fl) with the norm 

\\w\\bv = \\w\\ L i + TV(w). 

We refer the reader to [1] [28] for the most relevant facts about such functions. 

To carry out the finite element discretization, we introduce the number of points 
in each direction iVj £ N, i = 1, 2. The mesh size in each direction is 

hi = w i > l = 1 > 2 - 

The nodes of the mesh are 



z k .i = (khi,lh 2 ), fc = 0,iVi, l = 0,N 2 . 

The triangulation is 
where 

T Kl = {z £ R 2 : zlj < z 1 < zl +hl , zlj < z 2 < zj^ +1 } . 

Clearly Cl = UtetT- On the basis of such triangulation we define 

(2.2) X h = {w h £ C Q (Cl) : w h \ T £ Qi, VT e T} . 

Notice that, given the special structure of the mesh, if we denote by A^z the La- 
grange basis function associated with node Zk,i, then 

A M (z) = A^(z 1 )A 2 (z 2 ), 

where Ai, i — 1,2, denotes the one-dimensional Lagrange basis function associated 
with node z\. on a mesh of size hi. 
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3. TV DIMINISHING INTERPOLANT: PERIODIC DOMAINS 



In this Section we construct the TV diminishing interpolation operator in the 
torus, Q — S 1 x S , and prove some of its approximation properties. It is based 
on a weighted averaged Taylor polynomial as in 

HO HU but, to better exploit the 
symmetries related to our mesh, we use a very specific weighting function. 
To begin with, let xi be the characteristic function of the interval I and set 

= J^X[-hi/2,hi/2]> £ = 1,2. 
Notice that J ip l = 1 and HV^Hl 00 = 1/hi- The weighting function is then 

(3.i) ^oo = ^VtyV). 

Notice that 



hi 


hi 


X 


h 2 


hi 


2 


2 


2 


2 



(3.2) Q = supp-0 
Given a function w £ L 1 (f2), similarly to [71 [HI [12], we define 

(3.3) W ktl = f w(z Kl + C)^(C) d( G K. 



The interpolant HhW G Xh is the unique function that satisfies 
(3.4) n h w(z k>l ) = W kyl . 

The stability of this operator is as follows. 

Lemma 3.5 (instability). The operator Hh : L (CI) — > Xh is well defined and 
stable in any L p (fl), 1 < p < oo, namely 

(3.6) \\n h w\\ LP < \\w\\ LP , VweLP(n). 

Proof. It suffices to show (|3.6[) for p = 1 and p — oo . The intermediate powers will 
follow by function space interpolation. Estimate (|3.6[) for p = oo is evident. To 
obtain the ^-estimate we consider T K (£T and 



in^lU 1 ^) 



K+M+l 



J2 W k ,iAk,i(z) 



' k=K,l=l 

Notice that, if k = + 1 and I = £,£ + 1, 

/ A fe ./(z)dz= / 
By definition 



dz< \ w Ki\ / A fcl i(z)d«, 



A^dz 1 



\f(z 2 )dz 2 = ^hih 2 . 



\W k ,i\ < 



1 



HOIdC- 



hih 2 J ZkA+l _ 

From the considerations given above it follows that 

JVi-l,JV,-l , N 1 -l,N 2 -l K+1J+1 

E E 



|H h «;||£i = \\ u hw\\Li(T K , e ) < j 



HC)\dC 



K = 0^ = K=0,£=0 k = K,l=£ Jz k,l+Q 

To conclude it suffices to notice that every vertex belongs to exactly four cells. □ 
The main result of this contribution is the following. 
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Theorem 3.7 (TV diminishing: periodic case). The operator Hh : — > Xh 
defined in (|3.4j) satisfies (ITVDjl . 

Proof. Notice that, thanks to the particular structure of the mesh, the basis func- 
tions satisfy 

9xA M = (Ai)'A? = ±i-Af. 
If we C X (Q) and T = T M e T, then 

d 1 U h w\ T = Wu (^) Aj + Wfc+i,, (i) A 2 

(3-8) + W kJ+1 (^j (1 - A?) + W k+hl+1 (j^j (1 - A?) 

= 7^ {W k+U - W k ,i) A 2 + ^- (W k+U+1 - W fcl ,+i) (1 - A; 2 ), 

where we have used that, for every I, we have A 2 + Xf +1 = 1 on [z 2 v z\ In 
view of the definition of W k j, given in p. 31) . we can write 

W k +i,i-W ktl = f (w(z k+u +C) -w(z kA + ())^(C)dC 

JQ 

= hi / diw(sz k+u + (1 - s)z k ,i + CMC) dsd(. 
Notice that, at this point, it was essential to have a C 1 function. Notice also that 

\f(z 2 )dz = / {l-X\{z 2 ))dz=\\T\ = \h x h 2 . 
it Jt z z 

With these observations we estimate the local variation in the first coordinate di- 
rection of the interpolant as 



/ \d x n hW {z)\dz<\h x h 2 [ [ 

Jt z Jq Jo 



\d lW (sz k+hl + (1 - s)z k> i + C)| ^(C) dsdC 



+ ~h 1 h 2 [ [ \diw(sz k+lti+1 + (1 - s)z kA+1 + C)| V'(C) dsdC- 
z Jq Jo 

The total variation in the first coordinate direction is then 

/ \d 1 Il h w\= E / Iftnfc«;(2)|dz 
^ a ;=o fe=o Tfc - ; 

<~hih 2 V V / (|&u;(azfc + i,i + (l-*)«*,j+OI 

2 •/« i=o fc=0 Jo 
+ |#nu(sz fc+ i )i+1 + (1 - s)z fe! i + i + C)|) ^(0 ds dC- 

Introduce the change of variables (0, 1) 9 s i-» <7fc G ; + C\ ^+1 1 + C 1 ) given by 

(Tfe = + (! - sK.i + C\ 
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and notice that with it we can add over k to obtain 

/ \diU h w\<h 2 [ V ( I' \d lW (aX 2 + lh 2 )\da 
Jn 1 Jo ,_ n \Jo 



1=0 

1 

2 



+ / \d lW (a,C z + (l + l)h 2 )| daj^(C)dC 

^2~l r h 2 /2 



rhi/2 Na-1 r h 2 /2 r l 

<hJ w^E/ / 

J—hi/2 i=0 J-h 2 /2J0 



where, to arrive at the last inequality, we have used the definition of ip given in 
(f3~Tj) . Finally, we recall that ||^ 2 ||l°° = J "0 1 = 1, and 



W 2-l ph 2 /2 rl p 

V / / |9i W (a,C 2 +^2)|d ( 7dC 2 = / 

,_ n J-h 2 /2J0 Jn 



to conclude the proof for a C 1 function w. 

In general, i.e., w G .BV(f2) \C 1 (Q), we use an approximation argument. We 
recall [23 Theorem 5.3.3] that smooth functions are dense in BV(fl) under strict 
convergence. In other words, there is a sequence {w„} n gN C C°°(Q) such that 

lim \\w n — w\\li = 0, limsupTV(w„) < TV(iy). 

n ^oo n— >oo 

Since TV (HhW n ) < TV(w„) for all n E N, we deduce 

TV(n h w) = lim TV^Wn) < liminf TV(w„) < TV(u>), 

n— f oo n— > 

and thus conclude the proof. □ 
The approximation properties of this operator are summarized in the following. 
Theorem 3.9 (Approximation). Let h = ma,x{hi : i = 1, 2}. If w € £JV(r2), t/ien 

(3.10) ||w - n^wllii < ch TV(u>), 

wii/i c > a geometric constant. If w € _By(f2) n L p (f2) wzt/i p > 1, i/ien /or a/Z 

ge [l,p], 

(3.11) \\w - n fc «;|U, < ctf-'TVlwf-'WwWi,, i = i-i + -. 

Ifwe Wp(Ci) with 1 < p < oo, then 

(3.12) \\w - n h w|| LP < c/i 2 |w| W 2. 

Proof. For (|3.12p we use the symmetries of the mesh and the averaging procedure. 
In other words, if w G Pi, then w(zk,i) — Wk,i and HhW = w. This, together with 
an argument a la Bramble-Hilbert, implies the result. 

Estimate (|3.11[) follows from (13.61) and f|3 . 10[) using the well-known function space 
interpolation inequality [26] 

IMk« < IMIiilHIi,, - = i-^ + -. 

u q 1 p 
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We only sketch the proof of the first statement, as the arguments are similar to 
those used to obtain Theorem 13.71 We write 



Nt— l,iV2— l 



\\w - U h w\\ L i = E / \w-U h w\, 

and estimate the local differences as 

| to - U h w\ <\w- Wk a ,l \ + \W ko ,i - Uhw\, 

where (fccb'o) is such that Zfc 0l ; is a vertex of T Kj f. Provided that w £ C 1 (0) the 
bound on the first term is immediate, for if z € T K £, 



\w(z) - W ko ,i \ 



(w(z) - w(z ko j a + C)) ip(0 dC 



so that 



<ch / |V W ( S z + (l- S )(z feo , /o +C))|^(C)dsdC, 

JQ JO 

JVi-l,JV 2 -l 

E / \™-W ko , lo \<ch / |Vti;|, 



K,i=0 

because HV'IIl 00 = (^1^-2) 1 = I^k.^" 1 . For the second term we recall that the basis 
functions form a partition of unity, whence 



\W kotlo -U h w(z)\ 



W ko ,l Yl A kA z )- E W k,l A k,l( z ) 



k=K,l=l 



K+l,£+l 



< J2 A kM\w hl i-w ko ,i \ : 



Using, again, the differentiability properties of w, 

\W k .i-W ko .i \ <h f f |Vui (sz k .i + (1 - s)(z ka .i + C))| dsV(C) dC, 

JQ JO 

which, employing the facts that (A^fz)! < 1 and ||i/ , ||l° o = {hih2)~ 1 1 implies 

iVl-l,iV2-l 



E / \H h w-W ko<lo \<ch f \Vw\ 

. « r. JT K , JQ 



The proof concludes via a density argument with C 1 functions. 



□ 



4. TV DIMINISHING INTERPOLANT: NON-PERIODIC DOMAINS 

Let us describe two possible ways of modifying our construction in the case Q = 
(0, l) 2 . Both constructions will satisfy (jTVDj) and possess optimal approximation 
properties in the space BV(fl), so that all the results presented in Section [5] still 
hold in this case. Property (|3.12j) however, will no longer be valid in general. 
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4.1. Interpolation based on homotetic transformations of the domain. Set 

e = A max{hi, /12}, we define 

n e = (-e,l + ef 
and notice that fi <s f2 e . For a function w G £?V^(f2) we define 

( \ 1 ( v±L\ c o 

Notice that, employing the change of variables x = (y+e)/(l+2e) valid for functions 
of bounded variation [TJ Theorem 3.16], we obtain 

y + 



TVn.Ctfle* ' 



l + 2e 



(l + 2e) 2 
= TVn(io), 

where by TVa(w) we have denoted the total variation of the function w over the 
set A. Clearly w e \n € 51/(0). 

We define I^u; = II/ji^, where the operator II/,, is defined as in Section[3] Notice 
that the averaging procedure necessary for the definition of HhW e makes sense, since 
the function w e is defined in f2 £ . The results of Section [3] then imply 

TV Q (l h w) = TV Q (a h w e ) < TVnM = TV n (w). 

The approximation properties (|3.10|) and p.lip of 2/j can be obtained similarly. 

Remark 4.1 (Second order approximation). The approximation property (|3.12[) 
cannot hold for X/j. In fact if w = c, then I^w — c/(l + 2e). This also shows that 
the operator does not preserve constants, which might be an undesirable feature. 

4.2. Interpolation based on rescaled local averages. The construction of §4.1[ 
although simple, does not preserve constants. For this reason we consider a slightly 
more complicated procedure. It is a modification of the construction of Section [3] 
to take into account the effects of the boundary. 

The operator Ch '■ L 1 (il) 3 w 1— > ChW £ Xh is then such that if Zk,i is an internal 
node, the value of the interpolant is as in (j3.4|) . On the other hand, if Zk,i is a 
boundary node, we define 

Qn(M) = (z k j + Q)nn, 

and 

(4-2) W k<l = 1 / w(y)dy. 

\Qn{k,l)\ J Qn (kd) 

Observe that, \Q\/\Qn(k,l)\ = 4 for vertices of il and |Q|/|Qn(M)l = 2 f° r tne 
rest of boundary nodes. This rescaling allows us to take into account the lack of 
symmetry on the interaction with neighboring nodes. 

Theorem 4.3 (TV diminishing: non-periodic case). The operator Ch ■ L 1 (fi) — > 
Xh defined by Chw(zk,i) = Wk.i with W^.i as in (|3.3p for k — 1,N% — 1, I = 
1,N2 — 1 and by (|4.2|) otherwise, satisfies (jTVDp . 

Proof. We split Q into four congruent rectangles 

Q++ = {z e Q : z 1 > 0, z 2 > 0}, Q+- = {z e Q : z 1 > 0, z 2 < 0}, 
Q~+ = {z e Q : z 1 < 0, z 2 > 0}, Q— = {z e Q : z 1 < 0, z 2 < 0}, 
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and introduce 



Q~ + UQ H 



i? 4 



Q+- u Q H 



It is necessary to compute the contribution of the cells that intersect the bound- 
ary. To illustrate the procedure we consider Tq q and recall (|3.8p to write: 



2 "'Tb ) o 



\diC h w\ < 



1 



R- 



w(zi,o + C) dC - i^. 



1 



w\J Q w{zl < l+0dC -w-\L 



The first two terms on the right hand side can be estimated as 



hih 2 Jr-+ 



w(zx, + Q dC- 



h Y h 2 



u>(zo,o + C) dC 



2 

< — 

" h. 



2 JQ++ Jo 



2 JQ+ + 



|3iw(sz ,i + (1 - s)z ,o + 01 dCds 



2 JQ-+ JO 



\diw (szi.o + (1 - s)(^o,o + ^ei + C)) | dCds. 



The summation procedure is now as before. For this it suffices to notice that 

\Q ++ \ = \Q + -\ = \Q- + \ = \Q-\ = k\Ql D 



The operator Ch is LP -stable as in ([376]) . Its approximation properties are as in 
(|3.10[) and p. lip . For brevity we skip the proofs, as they repeat the arguments we 
have already presented. Property (|3.12|) however, will not be valid in general. 

Remark 4.4 (Second order approximation). Denote 

M e = {x G CI : dist(x, dfl) < e} . 

The proof of Theorem [4731 shows that if w G W^Ne) and \h < e, then 

\\w - C h w\\ L ^ (Afh) < ch\w\w^(Af 3h/2 )- 
Therefore, since \Afh\ < ch, we deduce 

\\w — Chw\\^i < c/i 2 . 



5. Applications 

In this Section we present two applications of the operators constructed in Sec- 
tions [3] and |4j The fact that the operators satisfy (|TVD[) and (13.10[) allows us to 
not only improve on existing results but also to have simpler arguments. 

5.1. Total Variation Minimization. Let a > and / G L°°(Q), Set 

(5.1) E(w)=TV(w) + ^\\w-ff L2 , 

and consider the problem of finding v G BV(fl) that minimizes this functional. 
This problem arises, for instance, in connection with image processing [22] . In [5] 
S. Bartels presents a first order finite element method for the approximation of 
minimizcrs of (|5.1[) as well as an iterative scheme for their computation, thereby 
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extending prior work of A. Chambolle and T. Pock [6] and B. Lucier and J. Wang 
[27] . The main convergence result of this work reads as follows: Let 

(5.2) v = argmin E(w), 

w<EBV(n) 

and 

(5.3) v h = axgminE(w h ). 

w h ex h 

If v G B f ^{L 2 (n)), with P G (0, 1], then 

g 

(5.4) ||v — Vh\\z a < c/i^+W. 

Let us, with the help our TV diminishing operator, improve on this result, which 
is at best ©(/i 1 / 4 ). 

Theorem 5.5 (Total variation minimization). Let v be as in (|5.2j) and as in 
([53)) . ,4ssMme f/iaf v G BV{Vt) n tten 

||v-« & || za < >/2/i 1/2 (||v|| L( » + ||/|| L< ») 1/2 TV(v) 1 / 2 . 

Proof. Owing to the strict convexity of E and the fact that v% is a minimizer over 
Xh, we have 

|||v - v h \\l a < E(v h ) - £?(v) < E{n h v) - E{v) 

< TV(IL lV ) - TV(v) + f (||n h v - /||| 2 - ||v - /||| 2 ) . 
Since satisfies ([TVDI) and (ET51) . 

||v-wfc||i a < (n /lV - v, n /lV + v - 2/) <2||n /l v-v|| L i(||v|| L oc. + ||/|| L oo) 

<2ft(||v|| i - + ||/|| i o.)TV(v), 
where, to obtain the last inequality, we used Theorem 13.91 □ 

Remark 5.6 (Reduced regularity). The conclusion of Theorem 15.51 is indeed an 
improvement over (|5.4[) . It is not difficult to show that BV(Vt) n L°°(Cl) ^ 
B^ 2 (L 2 (r2)); see [26j Lemma 38.1] for a proof and, in some sense, the converse 
statement. In this case f)5 -4[) yields a ©(/i 1 / 6 ) error estimate. Our results yields a 
rate which, given the smoothness of v, is optimal. 

5.2. TV flows. The equation 

(5-7) Ut=diV (^ul)' 

supplemented with suitable initial and boundary conditions is known as the TV 
flow, since it can be interpreted as the L 2 -subgradient flow of the functional TV. 
We refer the reader to [51 [TTJ [TH] for the analysis of this problem and to [H [131 [H] 
for its discretization. Concerning applications this problem, supplemented with 
a smooth strictly convex p-Laplacian term, arises, e.g., in the modeling of grain 
boundary motion [19] : facet formation and evolution [15] : and electromigration 
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The precise mathematical meaning of (|5.7|) is as a subgradient flow. In other 
words, given uo € L 2 (£l), we seek a function u : [0,T] — > L 2 (f2) fl £?y(f2) such that 

J (u t , u - w) + TV(u) - TVO) < 0, Vroe L 2 (ft) n BV{fl), 
I u| t=0 = u . 



(5.8) 



In prior work [3] , we provided an error analysis for discretizations of (|5.8[) as well 
as a convergent solution scheme for it. We considered the fully discrete scheme: 
Find {u|}f =0 C X h such that 



(5.9) 




k+1 - w h ) + TV« +1 ) - TV(w h ) < 0, Vw h e X h) 



T 



where At > is a time-step, K — [T/ At] and Th is some interpolation operator. 
Our main convergence result stated that if Uo € L°°(fi) n i?V^(S7) and TV(I/jUo) < 
c < oo uniformly in /i, then 

II u - ua||l-(l») < ||u -IfcUo|U> + c (At 1/2 + /* 1/6 ) , 

where u/j denotes the piecewise linear (in time) interpolant of the sequence {u k }. 
This result, as [20, 23 show, is optimal with respect to time, given the regularity 
of the initial data. 

Let us, with the help of the TV diminishing operator of Sections [3] and U improve 
the convergence rate in space. 

Theorem 5.10 (Convergence of TV flows). Let u solve (|5.8|) . Assume that the 
initial data u € L°°(0) H BV(fi). Let {u|}f =0 C X ft solve (|5T9]l . wit/i it£ = II^Uo. 

(5.11) || u - fi/,|U- ( £») < c (At 1 / 2 + h 1 ^) . 

Proof. The error of the time discretization is obtained with standard techniques, 
cf. [20j [23]. For this reason we concentrate on the error of space discretization 
and compare the solution of (|5.9p with the sequence {u k }^ =0 C L 2 (Q) n BV(Q), 
solution of the semidiscrete problem 

U A~T~ > uk+1 ~ w ) + TV ( ufc+1 )- TV H e L 2 (o)n w(f>), 

M° = Uq. 



Set w = in this inequality and wn = Hh,u k+1 in (|5.9[) and add them. Denoting 

||u fe+1 -n,V +1 || L2 



e = U — ttfj, we obtain 



|e fc+1 ||| 2 + ||e fc+1 - e fc || 2 2 < ||e fc || 2 2 + 2At 



u k+1 

a h 



At 



where we used (|TVDjl to cancel the TV terms. 

Setting Wh = ut in (|5.9I) is not difhcult to show that At^ 1 \\u k i +1 — w^|| £ 2 < 
cTV(uo) uniformly in h. 

The results of [4] show that if u G L°°{n) n BV(fl), then the solution to the 
semidiscrete flow {u k } C L°°(f2) n £?y(0). Consequently, it suffices to invoke 
Theorem 13.91 and add over k. □ 
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Remark 5.12 (Error in space). At this stage one might ask whether an optimal 
rate of convergence Oih 1 / 2 ) can be obtained for (|5.9[) . Our methods show that this 
is only possible under higher regularity on the solution. To realize this the reader 
can repeat our analysis for the heat equation (where the arguments are simpler) 
and convince himself that optimality is a consequence of the higher regularity that 
is expected from the solution. Unfortunately, the functional TV does not have any 
smoothing properties beyond the ones we have used so it not possible to obtain a 
better rate of convergence with the present techniques. 
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